Ultrafest
Stefano Coretta
03/04/2017
1 Data
num.splines <- 42
first.columns <- c(
"speaker",
"seconds",
"rec.date",
"prompt",
"label",
"TT.displacement.sm",
"TT.velocity",
"TT.velocity.abs",
"TD.displacement.sm",
"TD.velocity",
"TD.velocity.abs"
)
columns <- c(first.columns,
paste0(rep(c("X", "Y"), num.splines),
"_",
rep(1:num.splines, each = 2)
)
)
splines.raw <- list.files(path = "./voicing-effect/results",
pattern = "*-aaa.csv",
full.names = TRUE) %>%
map_df(~read.delim(., col.names = columns, na.strings = "*"))## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to charactersplines.closure.raw <- list.files(path = "./voicing-effect/results",
pattern = "*-closure.csv",
full.names = TRUE) %>%
map_df(~read.delim(., col.names = columns, na.strings = "*"))## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to charactersplines.phantom.raw <- list.files(path = "./voicing-effect/results",
pattern = "*-phantom.csv",
full.names = TRUE) %>%
map_df(~read_tsv(., col_names = columns, na = "*"))## Parsed with column specification:
## cols(
## .default = col_double(),
## speaker = col_character(),
## rec.date = col_character(),
## prompt = col_character(),
## label = col_character()
## )## See spec(...) for full column specifications.## Parsed with column specification:
## cols(
## .default = col_double(),
## speaker = col_character(),
## rec.date = col_character(),
## prompt = col_character(),
## label = col_character()
## )## See spec(...) for full column specifications.rm(num.splines, first.columns, columns)
languages <- read_csv("./voicing-effect/stimuli/languages.csv")## Parsed with column specification:
## cols(
## speaker = col_character(),
## language = col_character()
## )words <- read_csv("./voicing-effect/stimuli/nonce.csv")## Parsed with column specification:
## cols(
## item = col_integer(),
## word = col_character(),
## ipa = col_character(),
## c1 = col_character(),
## c1phonation = col_character(),
## vowel = col_character(),
## anteropost = col_character(),
## height = col_character(),
## c2 = col_character(),
## c2phonation = col_character(),
## c2place = col_character(),
## language = col_character()
## )splines <- splines.raw %>%
gather(spline, coordinate, matches("[XY]_")) %>%
separate(spline, c("axis", "fan"), convert = TRUE) %>%
spread(axis, coordinate) %>%
mutate(word = word(prompt, 2)) %>%
left_join(y = languages) %>%
left_join(y = words) %>%
mutate_if(is.character, as.factor)## Joining, by = "speaker"## Joining, by = c("word", "language")splines.it <- filter(splines, language == "italian")
splines.pl <- filter(splines, language == "polish")
splines.closure <- splines.closure.raw %>%
gather(spline, coordinate, matches("[XY]_")) %>%
separate(spline, c("axis", "fan"), convert = TRUE) %>%
spread(axis, coordinate) %>%
mutate(word = word(prompt, 2)) %>%
left_join(y = languages) %>%
left_join(y = words) %>%
mutate_if(is.character, as.factor)## Joining, by = "speaker"
## Joining, by = c("word", "language")splines.closure.it <- filter(splines.closure, language == "italian")
splines.closure.pl <- filter(splines.closure, language == "polish")2 Italian
2.1 Plot at maximum constriction
filter(splines.it, label == "max_TT") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel) +
scale_colour_discrete(name = "Phonation of C2")## Warning: Removed 754 rows containing non-finite values (stat_smooth).
filter(splines.it, label == "max_TD") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel) +
scale_colour_discrete(name = "Phonation of C2")## Warning: Removed 1127 rows containing non-finite values (stat_smooth).
filter(splines.it, label == "max_TT") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 754 rows containing non-finite values (stat_smooth).
filter(splines.it, label == "max_TD") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 1127 rows containing non-finite values (stat_smooth).
2.2 Plot at peak 1
filter(splines.it, label == "peak1_TT") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 844 rows containing non-finite values (stat_smooth).
filter(splines.it, label == "peak1_TD") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 997 rows containing non-finite values (stat_smooth).
filter(splines.it, label == "peak1_TT") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 844 rows containing non-finite values (stat_smooth).
filter(splines.it, label == "peak1_TD") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 997 rows containing non-finite values (stat_smooth).
2.3 Plot at closure
filter(splines.closure.it, c2place == "coronal") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 937 rows containing non-finite values (stat_smooth).
filter(splines.closure.it, c2place == "velar") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 1150 rows containing non-finite values (stat_smooth).
filter(splines.closure.it, c2place == "coronal") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 937 rows containing non-finite values (stat_smooth).
filter(splines.closure.it, c2place == "velar") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 1150 rows containing non-finite values (stat_smooth).
2.4 GAMMs at maximum constriction
splines.it$c2phonation.ord <- as.ordered(splines.it$c2phonation)
splines.it$c2place.ord <- as.ordered(splines.it$c2place)
splines.it$vowel.ord <- as.ordered(splines.it$vowel)
contrasts(splines.it$c2phonation.ord) <- "contr.treatment"
contrasts(splines.it$c2place.ord) <- "contr.treatment"
contrasts(splines.it$vowel.ord) <- "contr.treatment"
max.it <- filter(splines.it, label %in% c("max_TT", "max_TD"), c1phonation == "voiceless", vowel != "u") %>%
na.omit() %>%
arrange(rec.date, fan) %>%
start_event(column = "fan", event = c("rec.date"))max.it.gam <- bam(
Y ~
c2phonation.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr"),
data = max.it,
method = "ML"
)
summary(max.it.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.3889 0.1275 -42.259 < 2e-16 ***
## c2phonation.ord.L 1.2309 0.1809 6.804 1.25e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.636 8.938 923.974 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.817 3.523 9.283 2.33e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.849 Deviance explained = 85%
## -ML = 8962.5 Scale est. = 43.808 n = 2702max.it.gam.null <- bam(
Y ~
s(X, bs = "cr"),
data = max.it,
method = "ML"
)
compareML(max.it.gam, max.it.gam.null)## max.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## max.it.gam.null: Y ~ s(X, bs = "cr")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.it.gam.null 8999.026 3
## 2 max.it.gam 8962.513 6 36.514 3.000 9.592e-16 ***
##
## AIC difference: -70.19, model max.it.gam has lower AIC.plot_smooth(max.it.gam, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * X : numeric predictor; with 30 values ranging from -42.260800 to 59.923300.
plot_diff(max.it.gam, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * X : numeric predictor; with 100 values ranging from -42.260800 to 59.923300.
##
## X window(s) of significant difference(s):
## -42.260800 - -5.102945acf_plot(resid(max.it.gam), split_by=list(max.it$rec.date))
max.it.gam.2 <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = max.it,
method = "fREML"
)
summary(max.it.gam.2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.46938 0.08895 -61.490 <2e-16 ***
## c2phonation.ord.L 1.30798 0.12111 10.800 <2e-16 ***
## c2place.ord.L 1.13703 0.12529 9.075 <2e-16 ***
## vowel.ord.L 1.30514 0.12075 10.809 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.754 8.945 1013.76 <2e-16 ***
## s(X):c2phonation.ordvoiceless 6.292 7.238 18.51 <2e-16 ***
## s(X):c2place.ordvelar 8.701 8.946 314.19 <2e-16 ***
## s(X):vowel.ordo 6.661 7.835 28.53 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.933 Deviance explained = 93.4%
## fREML = 7909 Scale est. = 19.508 n = 2702max.it.gam.null <- bam(
Y ~
c2phonation.ord +
c2place.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr"),
data = max.it,
method = "fREML"
)
compareML(max.it.gam.2, max.it.gam.null)## max.it.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## max.it.gam.null: Y ~ c2phonation.ord + c2place.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr") + s(X, by = c2place.ord, bs = "cr")
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.it.gam.null 8059.126 9
## 2 max.it.gam.2 7908.981 12 150.145 3.000 < 2e-16 ***
##
## AIC difference: -315.03, model max.it.gam.2 has lower AIC.plot_smooth(max.it.gam.2, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -42.260800 to 59.923300.
plot_diff(max.it.gam.2, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -42.260800 to 59.923300.
##
## X window(s) of significant difference(s):
## -38.132149 - -12.328084acf_plot(resid(max.it.gam.2), split_by=list(max.it$rec.date))
max.it.gamm <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = max.it,
method = "fREML"
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(max.it.gamm)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X,
## by = vowel.ord, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.8210 4.4120 -2.226 0.02611 *
## c2phonation.ord.L 1.4489 0.4473 3.239 0.00121 **
## c2place.ord.L 0.8747 0.4548 1.923 0.05454 .
## vowel.ord.L 1.2193 0.4522 2.696 0.00706 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 6.161 6.763 6.500 1.21e-07 ***
## s(X):c2phonation.ordvoiceless 5.441 6.330 4.058 0.000408 ***
## s(X,rec.date) 234.516 412.000 4.536 < 2e-16 ***
## s(X,speaker) 6.481 8.000 32.197 < 2e-16 ***
## s(X):c2place.ordvelar 8.811 8.953 197.900 < 2e-16 ***
## s(X):vowel.ordo 6.656 7.497 11.144 1.54e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.965 Deviance explained = 96.9%
## fREML = 7334.8 Scale est. = 10.057 n = 2702max.it.gam.null <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = max.it,
method = "fREML"
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.compareML(max.it.gamm, max.it.gam.null)## max.it.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X,
## by = vowel.ord, bs = "cr")
##
## max.it.gam.null: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord)
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.it.gam.null 7430.089 15
## 2 max.it.gamm 7334.843 18 95.246 3.000 < 2e-16 ***
##
## AIC difference: -141.27, model max.it.gamm has lower AIC.plot_smooth(max.it.gamm, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -42.260800 to 59.923300.
## * rec.date : factor; set to the value(s): 29/11/2016 15:11:03.
## * speaker : factor; set to the value(s): sc01.
plot_diff(max.it.gamm, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -42.260800 to 59.923300.
## * rec.date : factor; set to the value(s): 29/11/2016 15:11:03.
## * speaker : factor; set to the value(s): sc01.
##
## X window(s) of significant difference(s):
## -42.260800 - -15.424572acf_plot(resid(max.it.gamm), split_by=list(max.it$rec.date))
2.4.1 AR1 model
r1 <- start_value_rho(max.it.gamm)
max.it.gam.AR <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = max.it,
method = "fREML",
rho = r1,
AR.start = max.it$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(max.it.gam.AR)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X,
## by = vowel.ord, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.4983 3.6253 -2.068 0.03871 *
## c2phonation.ord.L 1.4251 0.3367 4.233 2.39e-05 ***
## c2place.ord.L 1.0267 0.3385 3.033 0.00244 **
## vowel.ord.L 1.0714 0.3364 3.184 0.00147 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 7.254 7.812 7.852 4.16e-09 ***
## s(X):c2phonation.ordvoiceless 7.172 8.020 7.913 1.42e-10 ***
## s(X,rec.date) 134.559 412.000 0.807 < 2e-16 ***
## s(X,speaker) 5.706 8.000 25.634 < 2e-16 ***
## s(X):c2place.ordvelar 8.709 8.936 135.867 < 2e-16 ***
## s(X):vowel.ordo 5.644 6.654 8.406 1.41e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.953 Deviance explained = 95.6%
## fREML = 5823.8 Scale est. = 7.9222 n = 2702max.it.gam.AR.null <- bam(
Y ~
# c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
# s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = max.it,
method = "fREML",
rho = r1,
AR.start = max.it$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.compareML(max.it.gam.AR.null, max.it.gam.AR)## max.it.gam.AR.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date,
## bs = "fs", xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord, bs = "cr")
##
## max.it.gam.AR: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X,
## by = vowel.ord, bs = "cr")
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.it.gam.AR.null 5861.708 15
## 2 max.it.gam.AR 5823.848 18 37.860 3.000 2.539e-16 ***## Warning in compareML(max.it.gam.AR.null, max.it.gam.AR): AIC is not
## reliable, because an AR1 model is included (rho1 = 0.725285, rho2 =
## 0.725285).acf_plot(resid(max.it.gam.AR), split_by=list(max.it$rec.date))
acf_resid(max.it.gam.AR, split_pred = "AR.start")
plot_smooth(max.it.gam.AR, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -42.260800 to 59.923300.
## * rec.date : factor; set to the value(s): 29/11/2016 15:11:03.
## * speaker : factor; set to the value(s): sc01.
plot_diff(max.it.gam.AR, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -42.260800 to 59.923300.
## * rec.date : factor; set to the value(s): 29/11/2016 15:11:03.
## * speaker : factor; set to the value(s): sc01.
##
## X window(s) of significant difference(s):
## -41.228637 - -17.4888972.5 GAMMs at closure
splines.closure.it$c2phonation.ord <- as.ordered(splines.closure.it$c2phonation)
splines.closure.it$c2place.ord <- as.ordered(splines.closure.it$c2place)
splines.closure.it$vowel.ord <- as.ordered(splines.closure.it$vowel)
contrasts(splines.closure.it$c2phonation.ord) <- "contr.treatment"
contrasts(splines.closure.it$c2place.ord) <- "contr.treatment"
contrasts(splines.closure.it$vowel.ord) <- "contr.treatment"
closure.it <- splines.closure.it %>%
filter(c1phonation == "voiceless", vowel != "u") %>%
na.omit() %>%
arrange(rec.date, fan) %>%
start_event(column = "fan", event = c("rec.date"))closure.it.gam <- bam(
Y ~
c2phonation.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr"),
data = closure.it,
method = "ML"
)
summary(closure.it.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.9307 0.1124 -52.776 < 2e-16 ***
## c2phonation.ord.L 0.7339 0.1591 4.614 4.14e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.676 8.927 1037.842 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 6.166 7.122 9.889 2.09e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.86 Deviance explained = 86.1%
## -ML = 9003.8 Scale est. = 35.076 n = 2806closure.it.gam.null <- bam(
Y ~
s(X, bs = "cr"),
data = closure.it,
method = "ML"
)
compareML(closure.it.gam, closure.it.gam.null)## closure.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## closure.it.gam.null: Y ~ s(X, bs = "cr")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 closure.it.gam.null 9039.474 3
## 2 closure.it.gam 9003.765 6 35.709 3.000 2.120e-15 ***
##
## AIC difference: -81.83, model closure.it.gam has lower AIC.plot_smooth(closure.it.gam, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * X : numeric predictor; with 30 values ranging from -45.980000 to 64.280100.
plot_diff(closure.it.gam, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * X : numeric predictor; with 100 values ranging from -45.980000 to 64.280100.
##
## X window(s) of significant difference(s):
## -35.956355 - -19.250279
## -5.885418 - 10.820658
## 27.526733 - 38.664117acf_plot(resid(closure.it.gam), split_by=list(closure.it$rec.date))
closure.it.gam.2 <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = closure.it,
method = "fREML"
)
summary(closure.it.gam.2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.90863 0.08051 -73.390 < 2e-16 ***
## c2phonation.ord.L 0.94698 0.10990 8.617 < 2e-16 ***
## c2place.ord.L 0.73290 0.11296 6.488 1.03e-10 ***
## vowel.ord.L 1.29048 0.11020 11.710 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.749 8.927 1109.85 <2e-16 ***
## s(X):c2phonation.ordvoiceless 7.178 8.029 23.34 <2e-16 ***
## s(X):c2place.ordvelar 8.736 8.959 277.75 <2e-16 ***
## s(X):vowel.ordo 7.730 8.440 48.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.934 Deviance explained = 93.5%
## fREML = 7989.8 Scale est. = 16.566 n = 2806closure.it.gam.null <- bam(
Y ~
# c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
# s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = closure.it,
method = "fREML"
)
compareML(closure.it.gam.2, closure.it.gam.null)## closure.it.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord, bs = "cr")
##
## closure.it.gam.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord, bs = "cr")
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 closure.it.gam.null 8102.335 9
## 2 closure.it.gam.2 7989.848 12 112.487 3.000 < 2e-16 ***
##
## AIC difference: -242.02, model closure.it.gam.2 has lower AIC.plot_smooth(closure.it.gam.2, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -45.980000 to 64.280100.
plot_diff(closure.it.gam.2, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -45.980000 to 64.280100.
##
## X window(s) of significant difference(s):
## -45.980000 - -43.752523
## -35.956355 - -17.022802
## -2.544203 - 16.389349acf_plot(resid(closure.it.gam.2), split_by=list(closure.it$rec.date))
closure.it.gamm <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = closure.it,
method = "fREML"
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(closure.it.gamm)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X,
## by = vowel.ord, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.8567 4.2518 -2.083 0.03735 *
## c2phonation.ord.L 1.0214 0.4111 2.485 0.01303 *
## c2place.ord.L 0.6041 0.4129 1.463 0.14358
## vowel.ord.L 1.3767 0.4126 3.337 0.00086 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 6.835 7.419 9.060 1.89e-11 ***
## s(X):c2phonation.ordvoiceless 7.272 8.048 8.385 2.50e-11 ***
## s(X,rec.date) 249.926 402.000 5.429 < 2e-16 ***
## s(X,speaker) 5.928 8.000 33.672 < 2e-16 ***
## s(X):c2place.ordvelar 8.801 8.958 205.734 < 2e-16 ***
## s(X):vowel.ordo 8.593 8.884 23.072 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.973 Deviance explained = 97.5%
## fREML = 7121.2 Scale est. = 6.9052 n = 2806closure.it.gamm.null <- bam(
Y ~
# c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
# s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = closure.it,
method = "fREML"
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.compareML(closure.it.gamm, closure.it.gamm.null)## closure.it.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X,
## by = vowel.ord, bs = "cr")
##
## closure.it.gamm.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date,
## bs = "fs", xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord, bs = "cr")
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 closure.it.gamm.null 7157.603 15
## 2 closure.it.gamm 7121.171 18 36.432 3.000 1.040e-15 ***
##
## AIC difference: 31.42, model closure.it.gamm.null has lower AIC.AIC(closure.it.gamm, closure.it.gamm.null)## df AIC
## closure.it.gamm 295.5252 13668.49
## closure.it.gamm.null 311.7325 13637.07plot_smooth(closure.it.gamm, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -45.980000 to 64.280100.
## * rec.date : factor; set to the value(s): 12/12/2016 14:44:52.
## * speaker : factor; set to the value(s): sc01.
plot_diff(closure.it.gamm, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -45.980000 to 64.280100.
## * rec.date : factor; set to the value(s): 12/12/2016 14:44:52.
## * speaker : factor; set to the value(s): sc01.
##
## X window(s) of significant difference(s):
## -40.411308 - -18.136540acf_plot(resid(closure.it.gamm), split_by=list(closure.it$rec.date))
2.5.1 AR1 model
r1 <- start_value_rho(closure.it.gamm)
closure.it.gam.AR <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = closure.it,
method = "fREML",
rho = r1,
AR.start = closure.it$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(closure.it.gam.AR)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X,
## by = vowel.ord, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.1777 3.8108 -2.146 0.031970 *
## c2phonation.ord.L 1.0379 0.3333 3.114 0.001867 **
## c2place.ord.L 0.7008 0.3345 2.095 0.036297 *
## vowel.ord.L 1.2586 0.3341 3.768 0.000168 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 7.185 7.753 7.214 4.45e-09 ***
## s(X):c2phonation.ordvoiceless 6.983 7.901 8.287 5.12e-11 ***
## s(X,rec.date) 186.323 402.000 1.730 < 2e-16 ***
## s(X,speaker) 5.740 8.000 32.534 < 2e-16 ***
## s(X):c2place.ordvelar 8.730 8.946 144.770 < 2e-16 ***
## s(X):vowel.ordo 7.860 8.526 17.379 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.965 Deviance explained = 96.8%
## fREML = 5533.1 Scale est. = 5.2451 n = 2806closure.it.gam.AR.null <- bam(
Y ~
# c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
# s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = closure.it,
method = "fREML",
rho = r1,
AR.start = closure.it$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.compareML(closure.it.gam.AR, closure.it.gam.AR.null)## closure.it.gam.AR: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X,
## by = vowel.ord, bs = "cr")
##
## closure.it.gam.AR.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date,
## bs = "fs", xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord, bs = "cr")
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 closure.it.gam.AR.null 5566.155 15
## 2 closure.it.gam.AR 5533.075 18 33.079 3.000 2.835e-14 ***## Warning in compareML(closure.it.gam.AR, closure.it.gam.AR.null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.725722, rho2 =
## 0.725722).acf_plot(resid(closure.it.gam.AR), split_by=list(closure.it$rec.date))
acf_resid(closure.it.gam.AR, split_pred = "AR.start")
plot_smooth(closure.it.gam.AR, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -45.980000 to 64.280100.
## * rec.date : factor; set to the value(s): 12/12/2016 14:44:52.
## * speaker : factor; set to the value(s): sc01.
plot_diff(closure.it.gam.AR, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -45.980000 to 64.280100.
## * rec.date : factor; set to the value(s): 12/12/2016 14:44:52.
## * speaker : factor; set to the value(s): sc01.
##
## X window(s) of significant difference(s):
## -40.411308 - -15.9090642.6 GAMMs at peak 1
splines.it$c2phonation.ord <- as.ordered(splines.it$c2phonation)
contrasts(splines.it$c2phonation.ord) <- "contr.treatment"
peak.it <- filter(splines.it, label %in% c("peak1_TT", "peak1_TD"))
peak.it.gam <- bam(
Y ~
c2phonation.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr"),
data = peak.it,
method = "ML"
)
summary(peak.it.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.8024 0.1809 -32.083 <2e-16 ***
## c2phonation.ordvoiceless 0.5019 0.2541 1.976 0.0482 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.691 8.958 562.990 <2e-16 ***
## s(X):c2phonation.ordvoiceless 1.695 2.137 0.932 0.425
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.59 Deviance explained = 59%
## -ML = 22127 Scale est. = 96.198 n = 5971peak.it.gam.null <- bam(
Y ~
s(X, bs = "cr"),
data = peak.it,
method = "ML"
)
compareML(peak.it.gam, peak.it.gam.null)## peak.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## peak.it.gam.null: Y ~ s(X, bs = "cr")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 peak.it.gam.null 22129.00 3
## 2 peak.it.gam 22126.83 6 2.171 3.000 0.227
##
## AIC difference: -0.63, model peak.it.gam has lower AIC.## Warning in compareML(peak.it.gam, peak.it.gam.null): Only small difference in ML...plot_smooth(peak.it.gam, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * X : numeric predictor; with 30 values ranging from -48.091200 to 64.315100.
plot_diff(peak.it.gam, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * X : numeric predictor; with 100 values ranging from -48.091200 to 64.315100.
##
## Difference is not significant.2.7 Comparison of voiced consonants at closure and maximum constriction
voiced.it <- rbind(max.it, closure.it) %>%
filter(c2phonation == "voiced") %>%
mutate(label.ord = ifelse(label %in% c("max_TT", "max_TD"), "maximum", "closure"))
voiced.it$label.ord <- as.ordered(voiced.it$label.ord)
contrasts(voiced.it$label.ord) <- "contr.treatment"filter(voiced.it, c2place == "coronal") %>%
ggplot(aes(X, Y, colour = label.ord)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)
voiced.it.gamm <- bam(
Y ~
label.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = label.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = voiced.it,
method = "fREML"
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(voiced.it.gamm)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ label.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + s(X,
## by = label.ord, bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", m = 1,
## k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord,
## bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.3054 4.7990 -2.147 0.0319 *
## label.ordmaximum -1.3368 0.1183 -11.303 <2e-16 ***
## c2place.ord.L 0.9681 0.7984 1.212 0.2255
## vowel.ord.L 1.3332 0.7917 1.684 0.0923 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 5.881 6.501 7.462 1.32e-08 ***
## s(X):label.ordmaximum 7.188 8.018 69.108 < 2e-16 ***
## s(X,rec.date) 154.503 209.000 8.621 < 2e-16 ***
## s(X,speaker) 6.652 8.000 27.102 < 2e-16 ***
## s(X):c2place.ordvelar 8.742 8.929 90.637 < 2e-16 ***
## s(X):vowel.ordo 7.535 8.218 12.681 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.966 Deviance explained = 96.8%
## fREML = 7583.5 Scale est. = 9.5337 n = 2846r1 <- start_value_rho(voiced.it.gamm)
voiced.it.gam.AR <- bam(
Y ~
label.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = label.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = voiced.it,
method = "fREML",
rho = r1,
AR.start = voiced.it$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(voiced.it.gam.AR)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ label.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + s(X,
## by = label.ord, bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", m = 1,
## k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord,
## bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.4446 4.2142 -2.478 0.0133 *
## label.ordmaximum -1.2193 0.2554 -4.774 1.9e-06 ***
## c2place.ord.L 0.9768 0.4913 1.988 0.0469 *
## vowel.ord.L 1.0384 0.4902 2.118 0.0342 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 5.901 6.590 6.529 1.76e-07 ***
## s(X):label.ordmaximum 7.228 8.101 17.534 < 2e-16 ***
## s(X,rec.date) 109.376 209.000 2.344 < 2e-16 ***
## s(X,speaker) 6.672 8.000 25.111 < 2e-16 ***
## s(X):c2place.ordvelar 8.673 8.924 125.385 < 2e-16 ***
## s(X):vowel.ordo 7.029 7.916 13.625 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.957 Deviance explained = 95.9%
## fREML = 5914.4 Scale est. = 7.3839 n = 2846plot_smooth(voiced.it.gam.AR, view = "X",
plot_all= "label.ord", rug = FALSE)## Summary:
## * label.ord : factor; set to the value(s): closure, maximum.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -42.852800 to 64.280100.
## * rec.date : factor; set to the value(s): 29/11/2016 15:21:30.
## * speaker : factor; set to the value(s): sc01.
plot_diff(voiced.it.gam.AR, view = "X",
comp = list(label.ord = c("maximum", "closure")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -42.852800 to 64.280100.
## * rec.date : factor; set to the value(s): 29/11/2016 15:21:30.
## * speaker : factor; set to the value(s): sc01.
##
## X window(s) of significant difference(s):
## -38.524198 - -15.799037
## -2.813231 - 36.144187acf_resid(voiced.it.gam.AR, split_pred = "AR.start")
3 Polish
3.1 Plot at maximum constriction
filter(splines.pl, label == "max_TT") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1.5) +
facet_grid(speaker ~ vowel) +
scale_colour_discrete(name = "Phonation of C2")## Warning: Removed 274 rows containing non-finite values (stat_smooth).
filter(splines.pl, label == "max_TD") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel) +
scale_colour_discrete(name = "Phonation of C2")## Warning: Removed 328 rows containing non-finite values (stat_smooth).
filter(splines.pl, label == "max_TT") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1.5) +
facet_grid(speaker ~ vowel)## Warning: Removed 274 rows containing non-finite values (stat_smooth).
filter(splines.pl, label == "max_TD") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1.5) +
facet_grid(speaker ~ vowel)## Warning: Removed 328 rows containing non-finite values (stat_smooth).
3.2 Plot at peak 1
filter(splines.pl, label == "peak1_TT") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1.5) +
facet_grid(speaker ~ vowel)## Warning: Removed 273 rows containing non-finite values (stat_smooth).
filter(splines.pl, label == "peak1_TD") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1.5) +
facet_grid(speaker ~ vowel)## Warning: Removed 302 rows containing non-finite values (stat_smooth).
filter(splines.pl, label == "peak1_TT") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1.5) +
facet_grid(speaker ~ vowel)## Warning: Removed 273 rows containing non-finite values (stat_smooth).
filter(splines.pl, label == "peak1_TD") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1.5) +
facet_grid(speaker ~ vowel)## Warning: Removed 302 rows containing non-finite values (stat_smooth).
3.3 Plot at closure
filter(splines.closure.pl, c2place == "coronal") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1.5) +
facet_grid(speaker ~ vowel)## Warning: Removed 247 rows containing non-finite values (stat_smooth).
filter(splines.closure.pl, c2place == "velar") %>%
ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 284 rows containing non-finite values (stat_smooth).
filter(splines.closure.pl, c2place == "coronal") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 247 rows containing non-finite values (stat_smooth).
filter(splines.closure.pl, c2place == "velar") %>%
ggplot(aes(X, Y, colour = c2phonation)) +
geom_line(stat = "smooth", method = "loess") +
coord_fixed(ratio = 1) +
facet_grid(speaker ~ vowel)## Warning: Removed 284 rows containing non-finite values (stat_smooth).
3.4 GAMMs at maximum constriction ab03
splines.pl$c2phonation.ord <- as.ordered(splines.pl$c2phonation)
splines.pl$c2place.ord <- as.ordered(splines.pl$c2place)
splines.pl$vowel.ord <- as.ordered(splines.pl$vowel)
contrasts(splines.pl$c2phonation.ord) <- "contr.treatment"
contrasts(splines.pl$c2place.ord) <- "contr.treatment"
contrasts(splines.pl$vowel.ord) <- "contr.treatment"
max.pl <- splines.pl %>%
filter(label %in% c("max_TT", "max_TD"),
speaker == "ab03",
vowel != "u"
) %>%
na.omit() %>%
arrange(rec.date, fan) %>%
start_event(column = "fan", event = c("rec.date"))max.pl.gam <- bam(
Y ~
c2phonation.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr"),
data = max.pl,
method = "ML"
)
summary(max.pl.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9310 0.1236 7.530 8.09e-14 ***
## c2phonation.ord.L 0.7092 0.1752 4.048 5.40e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.865 8.993 332.184 <2e-16 ***
## s(X):c2phonation.ordvoiceless 2.129 2.662 3.826 0.0208 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.705 Deviance explained = 70.7%
## -ML = 5426.3 Scale est. = 26.88 n = 1763max.pl.gam.null <- bam(
Y ~
s(X, bs = "cr"),
data = max.pl,
method = "ML"
)
compareML(max.pl.gam, max.pl.gam.null)## max.pl.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## max.pl.gam.null: Y ~ s(X, bs = "cr")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.pl.gam.null 5437.967 3
## 2 max.pl.gam 5426.306 6 11.661 3.000 3.461e-05 ***
##
## AIC difference: -20.14, model max.pl.gam has lower AIC.plot_smooth(max.pl.gam, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * X : numeric predictor; with 30 values ranging from -34.964200 to 70.078500.
plot_diff(max.pl.gam, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500.
##
## X window(s) of significant difference(s):
## 23.392856 - 70.078500max.pl.gam.2 <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = max.pl,
method = "fREML"
)
summary(max.pl.gam.2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.86975 0.08813 9.869 < 2e-16 ***
## c2phonation.ord.L 0.63403 0.12287 5.160 2.75e-07 ***
## c2place.ord.L -0.17665 0.12468 -1.417 0.157
## vowel.ord.L -0.10397 0.12277 -0.847 0.397
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.660 8.929 204.37 <2e-16 ***
## s(X):c2phonation.ordvoiceless 2.511 3.135 3.37 0.0168 *
## s(X):c2place.ordvelar 8.412 8.879 174.24 <2e-16 ***
## s(X):vowel.ordo 8.275 8.848 32.24 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.856 Deviance explained = 85.8%
## fREML = 4825.9 Scale est. = 13.168 n = 1763max.pl.gam.null <- bam(
Y ~
# c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
# s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = max.pl,
method = "fREML"
)
compareML(max.pl.gam.2, max.pl.gam.null)## max.pl.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## max.pl.gam.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.pl.gam.null 4842.954 9
## 2 max.pl.gam.2 4825.858 12 17.095 3.000 1.806e-07 ***
##
## AIC difference: -31.01, model max.pl.gam.2 has lower AIC.plot_smooth(max.pl.gam.2, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -34.964200 to 70.078500.
plot_diff(max.pl.gam.2, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500.
##
## X window(s) of significant difference(s):
## -30.720051 - -11.621378
## 20.209743 - 70.078500max.pl.gamm <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
# s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = max.pl,
method = "fREML"
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(max.pl.gamm)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8847 0.5885 1.503 0.133
## c2phonation.ord.L 1.2735 0.7781 1.637 0.102
## c2place.ord.L -0.2137 0.8299 -0.257 0.797
## vowel.ord.L -0.1053 0.8298 -0.127 0.899
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.536 8.832 49.903 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 6.022 6.993 2.801 0.00728 **
## s(X,rec.date) 202.133 227.000 26.452 < 2e-16 ***
## s(X):c2place.ordvelar 8.749 8.933 30.231 < 2e-16 ***
## s(X):vowel.ordo 8.748 8.931 18.211 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.968 Deviance explained = 97.3%
## fREML = 3926.3 Scale est. = 2.8886 n = 1763max.pl.gamm.null <- bam(
Y ~
# c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
# s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
# s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = max.pl,
method = "fREML"
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.compareML(max.pl.gamm, max.pl.gam.null)## max.pl.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord, bs = "cr")
##
## max.pl.gam.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.pl.gam.null 4842.954 9
## 2 max.pl.gamm 3926.300 15 916.653 6.000 < 2e-16 ***
##
## AIC difference: -2517.12, model max.pl.gamm has lower AIC.plot_smooth(max.pl.gamm, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -34.964200 to 70.078500.
## * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.
plot_diff(max.pl.gamm, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500.
## * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.
##
## X window(s) of significant difference(s):
## -34.964200 - -22.2317523.4.1 AR1 model
r1 <- start_value_rho(max.pl.gamm)
max.pl.gam.AR <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = max.pl,
method = "fREML",
rho = r1,
AR.start = max.pl$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(max.pl.gam.AR)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8682 0.5156 1.684 0.0924 .
## c2phonation.ord.L 1.1093 0.6731 1.648 0.0995 .
## c2place.ord.L -0.1025 0.7273 -0.141 0.8879
## vowel.ord.L -0.1336 0.7276 -0.184 0.8543
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.456 8.824 43.295 <2e-16 ***
## s(X):c2phonation.ordvoiceless 4.942 5.993 1.600 0.124
## s(X,rec.date) 187.491 227.000 8.812 <2e-16 ***
## s(X):c2place.ordvelar 8.672 8.920 27.668 <2e-16 ***
## s(X):vowel.ordo 8.729 8.936 15.496 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.962 Deviance explained = 96.7%
## fREML = 3256.1 Scale est. = 2.4815 n = 1763max.pl.gam.AR.null <- bam(
Y ~
# c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
# s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = max.pl,
method = "fREML",
rho = r1,
AR.start = max.pl$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.compareML(max.pl.gam.AR.null, max.pl.gam.AR)## max.pl.gam.AR.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date,
## bs = "fs", xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord, bs = "cr")
##
## max.pl.gam.AR: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord, bs = "cr")
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.pl.gam.AR.null 3265.122 12
## 2 max.pl.gam.AR 3256.061 15 9.062 3.000 4.148e-04 ***## Warning in compareML(max.pl.gam.AR.null, max.pl.gam.AR): AIC is not
## reliable, because an AR1 model is included (rho1 = 0.615458, rho2 =
## 0.615458).acf_plot(resid(max.pl.gam.AR), split_by=list(max.pl$rec.date))
acf_resid(max.pl.gam.AR, split_pred = "AR.start")
plot_smooth(max.pl.gam.AR, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -34.964200 to 70.078500.
## * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.
plot_diff(max.pl.gam.AR, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500.
## * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.
##
## X window(s) of significant difference(s):
## 50.979827 - 54.1629393.5 GAMMs at maximum constriction ps02
max.pl.ps02 <- splines.pl %>%
filter(label %in% c("max_TT", "max_TD"),
speaker == "ps02",
vowel != "u",
X < 40,
X > -15
) %>%
na.omit() %>%
arrange(rec.date, fan) %>%
start_event(column = "fan", event = c("rec.date"))max.pl.gamm.ps02 <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
# s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = max.pl.ps02,
method = "fREML"
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(max.pl.gamm.ps02)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13.14685 0.18610 -70.644 < 2e-16 ***
## c2phonation.ord.L -0.08762 0.26224 -0.334 0.738325
## c2place.ord.L 1.83074 0.26315 6.957 4.8e-12 ***
## vowel.ord.L 0.96946 0.26082 3.717 0.000208 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.207 8.665 64.736 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 4.291 5.194 2.243 0.04978 *
## s(X,rec.date) 144.962 292.000 2.746 < 2e-16 ***
## s(X):c2place.ordvelar 8.804 8.951 90.910 < 2e-16 ***
## s(X):vowel.ordo 2.379 2.883 5.123 0.00166 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.863 Deviance explained = 87.5%
## fREML = 5726.5 Scale est. = 13.737 n = 20273.5.1 AR1 model
r1 <- start_value_rho(max.pl.gamm.ps02)
max.pl.gam.ps02.AR <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = max.pl.ps02,
method = "fREML",
rho = r1,
AR.start = max.pl.ps02$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(max.pl.gam.ps02.AR)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13.6822 0.2054 -66.597 < 2e-16 ***
## c2phonation.ord.L 0.1534 0.2793 0.549 0.582947
## c2place.ord.L 1.6141 0.2904 5.559 3.1e-08 ***
## vowel.ord.L 0.8489 0.2306 3.681 0.000239 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.163 8.667 65.792 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 4.255 5.311 1.418 0.219832
## s(X,rec.date) 125.209 292.000 1.364 < 2e-16 ***
## s(X):c2place.ordvelar 8.602 8.883 41.201 < 2e-16 ***
## s(X):vowel.ordo 1.000 1.001 11.903 0.000572 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.828 Deviance explained = 83.8%
## fREML = 4344.2 Scale est. = 7.8385 n = 2027max.pl.gam.ps02.AR.null <- bam(
Y ~
# c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
# s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord, bs = "cr"),
data = max.pl.ps02,
method = "fREML",
rho = r1,
AR.start = max.pl.ps02$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.compareML(max.pl.gam.ps02.AR.null, max.pl.gam.ps02.AR)## max.pl.gam.ps02.AR.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date,
## bs = "fs", xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord, bs = "cr")
##
## max.pl.gam.ps02.AR: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord, bs = "cr")
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.pl.gam.ps02.AR.null 4347.654 12
## 2 max.pl.gam.ps02.AR 4344.226 15 3.428 3.000 0.077## Warning in compareML(max.pl.gam.ps02.AR.null, max.pl.gam.ps02.AR): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.733987, rho2 =
## 0.733987).## Warning in compareML(max.pl.gam.ps02.AR.null, max.pl.gam.ps02.AR): Only small difference in fREML...acf_plot(resid(max.pl.gam.ps02.AR), split_by=list(max.pl.ps02$rec.date))
acf_resid(max.pl.gam.ps02.AR, split_pred = "AR.start")
plot_smooth(max.pl.gam.ps02.AR, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): velar.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -14.986000 to 39.996900.
## * rec.date : factor; set to the value(s): 07/02/2017 16:30:56.
plot_diff(max.pl.gam.ps02.AR, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): velar.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -14.986000 to 39.996900.
## * rec.date : factor; set to the value(s): 07/02/2017 16:30:56.
##
## Difference is not significant.3.6 GAMMs at closure
splines.closure.pl$c2phonation.ord <- as.ordered(splines.closure.pl$c2phonation)
splines.closure.pl$c2place.ord <- as.ordered(splines.closure.pl$c2place)
splines.closure.pl$vowel.ord <- as.ordered(splines.closure.pl$vowel)
contrasts(splines.closure.pl$c2phonation.ord) <- "contr.treatment"
contrasts(splines.closure.pl$c2place.ord) <- "contr.treatment"
contrasts(splines.closure.pl$vowel.ord) <- "contr.treatment"
closure.pl <- splines.closure.pl %>%
filter(speaker == "ab03", vowel != "u", X > -15) %>%
na.omit() %>%
arrange(rec.date, fan) %>%
start_event(column = "fan", event = c("rec.date"))closure.pl.gam <- bam(
Y ~
c2phonation.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr"),
data = closure.pl,
method = "ML"
)
summary(closure.pl.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.24355 0.14312 22.662 <2e-16 ***
## c2phonation.ord.L 0.05709 0.20260 0.282 0.778
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 7.429 8.356 128.934 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 1.894 2.367 7.644 0.000309 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.523 Deviance explained = 52.7%
## -ML = 4461.8 Scale est. = 29.302 n = 1432closure.pl.gam.null <- bam(
Y ~
s(X, bs = "cr"),
data = closure.pl,
method = "ML"
)
compareML(closure.pl.gam, closure.pl.gam.null)## closure.pl.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## closure.pl.gam.null: Y ~ s(X, bs = "cr")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 closure.pl.gam.null 4469.756 3
## 2 closure.pl.gam 4461.755 6 8.000 3.000 0.001 **
##
## AIC difference: -11.18, model closure.pl.gam has lower AIC.plot_smooth(closure.pl.gam, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * X : numeric predictor; with 30 values ranging from -14.986800 to 71.049000.
plot_diff(closure.pl.gam, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * X : numeric predictor; with 100 values ranging from -14.986800 to 71.049000.
##
## X window(s) of significant difference(s):
## -7.165364 - 11.084655
## 40.632303 - 71.049000closure.pl.gamm <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
# s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = closure.pl,
method = "fREML"
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(closure.pl.gamm)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.12889 0.40861 7.657 3.91e-14 ***
## c2phonation.ord.L -0.00129 0.54255 -0.002 0.9981
## c2place.ord.L -1.03057 0.57581 -1.790 0.0737 .
## vowel.ord.L -0.72943 0.55364 -1.318 0.1879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.785 8.915 47.967 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 5.262 6.307 1.741 0.21884
## s(X,rec.date) 207.831 227.000 88.103 < 2e-16 ***
## s(X):c2place.ordvelar 8.647 8.899 34.922 < 2e-16 ***
## s(X):vowel.ordo 6.395 7.429 3.651 0.00219 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.987 Deviance explained = 98.9%
## fREML = 2424.9 Scale est. = 0.82697 n = 1432closure.pl.gamm.null <- bam(
Y ~
# c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
# s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
# s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = closure.pl,
method = "fREML"
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.compareML(closure.pl.gamm, closure.pl.gamm.null)## closure.pl.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord)
##
## closure.pl.gamm.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date,
## bs = "fs", xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 closure.pl.gamm.null 2430.334 12
## 2 closure.pl.gamm 2424.913 15 5.421 3.000 0.013 *
##
## AIC difference: -10.11, model closure.pl.gamm has lower AIC.AIC(closure.pl.gamm, closure.pl.gamm.null)## df AIC
## closure.pl.gamm 243.6748 4015.338
## closure.pl.gamm.null 240.2176 4025.451plot_smooth(closure.pl.gamm, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -14.986800 to 71.049000.
## * rec.date : factor; set to the value(s): 24/03/2017 14:49:08.
plot_diff(closure.pl.gamm, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -14.986800 to 71.049000.
## * rec.date : factor; set to the value(s): 24/03/2017 14:49:08.
##
## Difference is not significant.acf_plot(resid(closure.pl.gamm), split_by=list(closure.pl$rec.date))
3.6.1 AR1 model
r1 <- start_value_rho(closure.pl.gamm)
closure.pl.gam.AR <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
# s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = closure.pl,
method = "fREML",
rho = r1,
AR.start = closure.pl$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.summary(closure.pl.gam.AR)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.1253 0.4017 7.779 1.57e-14 ***
## c2phonation.ord.L -0.0335 0.5372 -0.062 0.9503
## c2place.ord.L -1.0248 0.5664 -1.809 0.0706 .
## vowel.ord.L -0.7378 0.5462 -1.351 0.1770
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.751 8.908 42.801 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 5.019 6.109 1.417 0.26593
## s(X,rec.date) 201.924 227.000 41.057 < 2e-16 ***
## s(X):c2place.ordvelar 8.600 8.896 31.216 < 2e-16 ***
## s(X):vowel.ordo 6.057 7.195 3.311 0.00231 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.986 Deviance explained = 98.8%
## fREML = 2135 Scale est. = 0.7632 n = 1432closure.pl.gam.AR.null <- bam(
Y ~
# c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
# s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
# s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = closure.pl,
method = "fREML",
rho = r1,
AR.start = closure.pl$start.event
)## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.compareML(closure.pl.gam.AR, closure.pl.gam.AR.null)## closure.pl.gam.AR: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord)
##
## closure.pl.gam.AR.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date,
## bs = "fs", xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 closure.pl.gam.AR.null 2139.752 12
## 2 closure.pl.gam.AR 2134.956 15 4.796 3.000 0.022 *## Warning in compareML(closure.pl.gam.AR, closure.pl.gam.AR.null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.446313, rho2 =
## 0.446313).## Warning in compareML(closure.pl.gam.AR, closure.pl.gam.AR.null): Only small difference in fREML...acf_plot(resid(closure.pl.gam.AR), split_by=list(closure.pl$rec.date))
acf_resid(closure.pl.gam.AR, split_pred = "AR.start")
plot_smooth(closure.pl.gam.AR, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -14.986800 to 71.049000.
## * rec.date : factor; set to the value(s): 24/03/2017 14:49:08.
plot_diff(closure.pl.gam.AR, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -14.986800 to 71.049000.
## * rec.date : factor; set to the value(s): 24/03/2017 14:49:08.
##
## Difference is not significant.